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SECTION  I 


INTRODUCTION 

The  present  Author  has  developed  an  Alternating  Direction  Implicit  (ADI) 
technique  for  the  calculation  of  viscous,  incompressible,  steady  flows  pa3t 
an  arbitrary,  two-dimensional  body^.  Such  an  approach  used  the  vorticity- 
stream  function  Navier-Stokes  equations  in  a  system  of  general  body-fitted 
coordinates.  The  governing  equations  were  parabolized  by  adding  a  relaxation¬ 
like  time  derivative  to  the  stream  function  equation,  linearized  in  time  and 

2 

solved  by  means  of  the  ADI  procedure  of  Douglas  and  Gunn  .  The  method  used 

second-order-accurate  finite  differences  and,  at  convergence,  provided  a 

second-order-accurate  approximation  to  the  steady  flow  of  interest.  The 

major  limitation  of  the  proposed  ADI  approach  was  due  to  its  use  of  central 

differences,  which  limited  its  applicability  to  low  Reynolds  number  flows,  or 

to  separation-free  high  Reynolds  number  flows.  A  first-order-accurate  method, 

using  windward  differences  for  the  convective  terms  in  the  equations,  although 

feasible  in  principle,  is  not  recommended,  insofar  as  the  effective  Reynolds 

number  of  the  numerical  solution  is  lowered  by  the  numerical  viscosity 

introduced  by  the  first-order-accurate  windward  differences.  A  more  stable, 

viscosity-free  numerical  technique  is  obtainable  by  using  windward  differences 

for  the  convective  terms,  which  are  evaluated  implicitly,  and  correcting  them 

3 

to  second-order-accurate  central  differences,  explicitly  ,  or,  more  simply, 

4 

by  employing  the  incremental  (delta)  form  of  the  equations  and  using  wind¬ 
ward  differences  for  the  incremental  variables  and  central  differences  for 
the  nonincremental  variables.  Both  approaches,  which  are  shown  here  to 
coincide,  have  the  desirable  property  of  lowering  the  effective  Reynolds 
number  of  the  pseudo-transient  problem  by  adding  numerical  viscosity,  which 


is,  however,  completely  removed  from  the  solution  as  convergence  is  achieved. 
The  second  approach,  using  the  incremental  variables,  is  employed  in  this 
paper  to  provide  an  improved  version  of  the  ADI  technique  given  in  Ref.  1, 
which  is  more  stable  and  capable  of  resolving  high  Reynolds  number  separated 
flows,  while  maintaining  the  second-order  accuracy  of  the  numerical-viscosity- 
free  central  differences,  at  convergence. 

Furthermore,  following  the  idea  of  Rubin  and  Khosla^,  it  is  possible 
and  very  straight-forward  indeed  to  obtain  a  fourth-order-accurate  spline 
ADI  method  by  means  of  a  spline  deferred-corrector  approach"*.  The  present 
paper  also  provides  a  simplified  spline  ADI  technique  for  the  vorticity- 
stream  function  Navier-Stokes  equations,  which  has  all  the  features  of  the 
aforementioned  improved  ADI  method.  In  particular,  incremental  variables 
are  used  and,  in  order  to  enhance  the  stability  of  the  method,  the  convective 
terms  are  approximated  by  first-order-accurate  windward  differences.  At  the 
end  of  each  two-sweep  ADI  cycle,  the  right  hand  side  (RHS)  of  the  difference 
equations,  which  is  already  second-order-accurate,  is  corrected  by  means  of 
a  spline  interpolating  procedure,  explicitly,  so  that,  at  convergence,  a 
fourth-order-accurate  approximation  to  the  steady  flow  of  interest  is  obtained. 

The  present  report  develops  as  follows:  In  Section  II  the  basic  ideas 

of  obtaining  second-  or  fourth-order  accuracy  at  convergence,  while  using 

first-order-accurate  windward  differences  to  stabilize  the  transient  pheno- 

3  3 

menon,  are  provided.  The  approach  of  Rubin  and  Khosla  ’  is  briefly  reviewed 
and  the  equivalent  simpler  and  more  elegant  approach  used  in  this  study  is 
presented.  In  Section  III  those  ideas  are  applied  to  the  ADI  numerical 
technique  previously  developed  by  the  Author^,  to  provide  an  improved,  more 
stable,  second-order-accurate  ADI  method  as  well  as  a  fourth-order-accurate 
spline  ADI  procedure. 

2 


Finally,  the  results  obtained  by  applying  the  present  techniques  to 
two  model  problems  (viscous  flow  between  two  concentric  circles  and  the 


classical  driven  cavity  flow)  as  well  as  to  a  problem  of  practical  interest 
(viscous  flow  in  a  channel  of  complex  geometry)  are  presented  in  Section  IV. 


SECTION  II 


THE  DEFERRED  CORRECTOR  APPROACH  OF  RUBIN  AND  KHOSLA 

Rubin  and  Khosla  have  presented  their  simplified  spline  technique^, 

as  applied  to  the  numerical  solution  of  the  steady  state  Burgers  equation, 

starting  from  the  unsteady  equation: 

u  +  c  u  =  \)u 
t  X  XX 

The  numerical  procedure  of  Rubin  and  Khosla  uses  the  following  discrete 
form  of  Eqn.  (1) : 


n+l  n  rt+1  n+l.,  ,,  „  ,  n+1 

u.  -  u,  +  cA(u.  -  u.  ,)k/h  -  B  vfu.,, 
3  j  J  3-1  3+1 


„  n+l  n+l..  .2 

2u.  +  u.  , )k/h  = 

3  3~1 


c  {A(un  -  un  .)  -  hmn}  k/h  +  V  k{Mn  -  B(u"  -  2  Un  +  un  ,)/h2}  (2) 

3  3~1  3  3  3+1  3  J-l 


where  the  subscript  j,  j-l  and  j+1  indicate  the  spatial  grid  locations,  the 
superscripts  n  and  n+l  indicate  the  old  and  new  time  levels,  k  is  the  time 
step,  h  is  the  spatial  meshwidth,  m  and  M  are  the  fourth-order-accurate 
spline  approximations  of  and  u  and  A  and  B  are  two  arbitrary  constants, 
necessary  for  the  stability  of  the  method,  see  Ref.  5  for  details.  Notice 
that  the  high-order  spline  correction  terms  only  appear  on  the  old  time 
level,  known  RHS  of  Eqn.  (2),  so  that  at  each  time  advancement  only  a  tri¬ 
diagonal  system  has  to  be  solved,  exactly  as  for  the  case  of  a  low-order- 
accurate  finite  difference  scheme,  but  that,  . t  convergence,  a  fourth-order- 
accurate  solution  is  obtained.  Also  notice  that,  if  m  and  M  are  replaced  with 
simple  second-order-accurate  finite  difference  approximations,  Eqn.  (2)  leads 
to  the  unconditionally  stable  KR  scheme  (A  =  B  =  1)  of  Ref.  3. 

In  the  present  paper  an  ADI  and  a  spline  ADI  technique  will  be  developed 
for  the  Navier-Stokes  equations,  which  are  based  on  these  two  techniques. 


|i  ».i  i  in  i  ' 


u 


However,  an  Incremental  (delta)  formulation  is  preferred  In  this  paper  for 
its  superior  simplicity  and  elegance.  For  example,  by  writing  Eqn.  (2)  in 
delta  form,  one  obtains: 


where  Du  *=  un+^  -  u° .  Eqn.  (3)  actually  coincides  with  Eqn.  (2)  and  will 
obviously  produce  the  same  numerical  results,  but  it  is  simpler,  more  elegant 
and  requires  less  computation  effort.  Notice,  for  example,  that  the  RHS  of 
Eqn.  (3)  is  a  fourth-order-accurate  discrete  approximation  of  the  steady 
Burgers  equation,  which  thus,  at  convergence  (Du  =  0  at  all  gridpoints) ,  will 
be  resolved  with  the  desired  level  of  accuracy.  Moreover,  the  arbitrary 
constants  A  and  B,  which  only  multiply  the  incremental  variables,  are  clearly 
seen  not  to  influence  the  final  steady  state  solution.  By  comparing  the 
simplicity  of  Eqn.  (3)  with  respect  to  Eqn.  (2),  it  is  easy  to  understand 
the  advantage  of  using  the  delta  approach  in  complex  numerical  techniques 
for  the  Navier-Stokes  equations.  Obviously,  also  in  Eqn.  (3),  by  replacing 
m  and  M  with  standard  second-order-accurate  finite  differences,  one  obtains  a 
scheme  which  has  the  stability  of  a  windward  difference  scheme  and  the 
accuracy  (at  convergence)  of  a  central  difference  scheme.  In  such  a  case, 
although  unnecessary  for  stability,  the  use  of  values  greater  than  one  for 
A  and  B  can  further  enhance  the  convergence  of  the  numerical  method.  It 
remains  to  be  said  how  m  and  M  are  evaluated.  After  all  n+1  values  have 
been  obtained  directly  from  Eqn.  (2),  or  from  Eqn.  (3)  and  the  definition  of 


Du,  the  new  time  level  m  and  M  values  are  explicitly  evaluated  as: 


*  =  <uj+1  -  «J_1)/2  h  +  h(K  -  Kj+1)/12 


(4) 


and 


M.  =  {K.  +  (u...  -  2u.  +  u.  .)/h  }/2 
J  J  j+i  j  J-l 


(5) 


where  all  terms  are  easily  obtained  by  solving  the  following  tridiagonal 
system 


Vi +  4KJ +  Kj-1  ■  6<Vi  -  2UJ +  Vi)/h 


(6) 


see  Ref.  5  or  Ref.  8  for  details.  It  is  noteworthy  that,  with  respect  to  a 
standard  second-order-accurate  finite  difference  method,  the  simplified 
spline  technique  requires  the  solution  of  an  additional  tridiagonal  system 

g 

(Eqn.  6),  at  each  time  level,  whereas  a  standard  spline  technique  would 
require  the  solution  of  a  2  x  2  block-tridiagonal  system,  at  each  time 
level.  The  convenience  of  the  simplified  approach  is  seen  to  increase  when 
dealing  with  coupled  systems  of  equations  (e.g.  with  the  Navier-Stokes 
equations) . 


i 


t 
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SECTION  III 


THE  PRESENT  ADI  AND  SIMPLIFIED  SPLINE  ADI  METHODS 
FOR  THE  NAVIER-STOKES  EQUATIONS 

The  vorticity-stream  function  Navier-Stokes  equations  in  a  general 
system  of  curvilinear  body-oriented  coordinates  (£,n)  are  given ^  as: 

2 

<ot  +  (i/^oj^  -  J  -  (au^  -  2^n  +  Y“nT1  +  au^  +  tw^/J  Re  -  0  (7) 

and 

(«l»££  “  +  o^n  +  t^)/J2  +  u  -  t|»t  (8) 

where  a,  0,  y,  a  and  T  are  the  metric  coefficients  and  J  is  the  Jacobian  of 

the  coordinate  transformation,  and  a  relaxation-like  time  derivative  has 

been  added  to  the  stream  function  equation,  in  order  to  deal  with  a  parabolic 

system  of  equations^’ ^ .  For  this  reason,  Eqns.  (7)  and  (8)  are  not  the  time- 

dependent  Navier-Stokes  equations  and  in  the  present  time-marching  numerical 

techniques  only  the  converged  solutions  will  have  physical  meaning. 

Equations  (7,  8)  are  written  in  terms  of  the  incremental  variables, 

DtJj  =  -  (j/1.  Dm  =  wn+1  -  oin,  and  linearized  in  time  by  a  Taylor's  series 

2 

expansion,  which  neglects  terms  of  order  D  ,  to  give: 

Du/k  +  ( Dip  +  4>n  DUr  ~  Du)  )/J  -  (a  Du )frf.  +  y  Du  + 

a  Dun  +  T  Du)^)/J2Re  -  SSVEn  (9) 

and 

A 

D4»/k  -  Du  -  (a  Dip^  +  Y  D^rin  +  a  +  x  D^)/J  -  SSSFEn  (10) 

where  SSVEn  and  SSSFEn  are  shorthand  notations  for  steady  state  vorticity 
(stream  function)  equation  evaluated  at  the  n  time  level.  It  is  noteworthy 
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that  the  linearized  Eqns.  (9)  and  (10)  are  fully  implicit  for  the  incremental 
variables  except  for  the  mixed  derivatives,  which  are  treated  explicitly. 

Also,  for  more  generality,  two  constants,  A  and  B,  can  be  introduced  to 
multiply  all  convective  and  diffusive  incremental  terms,  respectively.  At 
this  point,  we  are  ready  to  solve  Eqns.  (9)  and  (10)  by  means  of  two-sweep 

ADI  techniques,  differing  only  in  the  level  of  accuracy  used  to  approximate 

SSVEn  and  SSSFEn.  These  ADI  methods,  derived  from  that  of  Douglas  and  Gunn2, 
proceed  as  follows:  at  the  first  sweep,  the  n  derivatives  and  the  source-like 
terms  in  the  LHS  (left  hand  side)  of  Eqns.  (9)  and  (10)  are  evaluated  implicitly 
to  give: 

Dw/k  +  (D^  oi-  -  l|/-Dw  )/J  -  (Y  Doj  +  a  Du  )/J2  Re  =  SSVEn  (11) 

n  S  K  n  nn  n 

Dpk  -  Du  -  (yD$  +  O  D^  )/J2  =  SSSFEn  (12) 

nn  n 

where  the  n,  indicate  that  the  solution  is  a  first  sweep  (predictor-type)  one 
and  all  the  £  derivatives,  which  are  evaluated  explicitly,  give  zero  contribu¬ 
tion  in  the  incremental  variables.  At  the  second  and  final  sweep,  all  the  £ 
derivatives  and  the  source-like  terms  are  evaluated  implicitly,  whereas  the 
n  derivatives  are  evaluated  at  the  first  sweep  (y)  level,  explicitly.  The 
resulting  equations  are  not  given  here  because,  for  computational  convenience, 
they  are  replaced  by  the  following  ones,  obtained  by  subtracting  from  ti»o*n 
the  first  sweep  Eqns.  (11,12): 

Dco/k  +  (t|/n  Du.  -  Dip.  un)/J  -  (a  Dow  +  T  Du).)/J2Re  =  Du/k  (13) 

n  s  s  n 

Dip/k  -  Du  -  (a  +  0  Dipg)/J2  *■  Dt^/k  -  Du  (14) 
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In  Eqns.  (11)  thru  (14)  all  the  incremental  terms  are  approximated  with 
central  differences  (the  second  derivatives)  and  windward  differences  (the 
first  derivatives)  whereas  the  nonincremental  terms  are  approximated  with 
standard  central  differences  (for  the  case  of  the  ADI  method)  or  fourth- 
order-accurate  spline  approximations  (for  the  case  of  the  spline  ADI  method). 
Therefore,  for  both  techniques,  a  series  of  2  x  2  block-tridiagonal  systems 
is  to  be  solved  at  each  sweep  of  the  ADI  procedure,  exactly  as  in  the  former 
ADI  technique  due  to  the  Author''". 

At  the  end  of  a  complete  ADI  cycle  the  solution  is  updated  as: 

i//n+1  =  <pn  +  Di//  (15) 

u)n+1  =  o)n  +  Do:  (16) 

and,  for  the  case  of  the  ADI  method,  the  process  is  repeated  until  a  satis¬ 
factory  convergence  is  achieved. 

For  the  case  of  the  spline  ADI,  however,  in  order  to  be  able  to  evaluate 
the  RHS  in  Eqns.  (11)  and  (12)  with  fourth-order  accuracy  (at  convergence), 
it  is  necessary  to  obtain  a  fifth-order  interpolating  polynomial  approximat¬ 
ing  the  new  values  of  the  stream  function  and  of  the  vorticity  along  each  row 
and  column  of  the  computational  grid.  This  is  done  by  solving  two  tridiagonal 
systems,  formally  identical  to  Eqn.  (6),  for  each  row  and  for  each  column  of 

gridpoints,  which  allow  to  evaluate  the  four  matrices  ,  Ko)£.  . 

1  >3  1 » J  1 » 3 

and  Kwn,  from  these,  all  the  corresponding  first  and  second  derivatives, 

1  *  J 

mi p£.  ,,  M ip£  ,,  etc.,  can  then  be  evaluated  explicitly,  by  means  of  expressions 
*  *  J  J 

formally  identical  to  Eqns.  (4)  and  (5).  As  far  as  the  boundary  conditions  are 
concerned,  for  the  case  of  Eqn.  (6),  the  first  and  last  values  of  K  are 


evaluated  by  linear  or  quadratic  extrapolations  from  the  neighboring  points. 

It  is  noteworthy  that  in  the  present  spline  ADI  technique,  if  a  N  x  N  mesh 
is  used,  at  each  ADI  sweep  it  is  necessary  to  solve  2N  2  x  2  block-tridiagonal 
systems  (Eqns.  7-10)  and  4N  simple  tridiagonal  systems  (equations  of  the  type 
of  Eqn.  6).  A  standard  spline  ADI  would  require  the  solution  of  2N  4  x  4 
block-tridiagonal  systems,  i.e.,  a  lot  more  computational  work.  Also,  the 
additional  work  of  the  present  simplified  spline  ADI  method  with  respect  to 
the  corresponding  ADI  approach  is  minimal,  considering  the  accuracy  improve¬ 
ment  it  provides. 

The  present  approaches  solve,  at  each  sweep,  the  vorticity  and  the 
stream  function  equations  as  a  coupled  set  on  each  row  and  column  of  the 
computational  grid.  For  this  reason  it  is  possible,  in  the  solution  routine 
for  each  2x2  block-tridiagonal  system,  to  accommodate  the  double  specifica¬ 
tion  on  the  stream  function  at  the  boundary  and  to  evaluate  the  vorticity 
at  the  wall,  directly.  For  the  case  of  the  ADI  method,  the  boundary  condi¬ 
tions  are  imposed  exactly  as  in  Ref.  1,  with  the  difference  that,  in  the 
present  case,  the  incremental  approach  is  used  also  at  the  first  sweep  of  the 
ADI  method.  For  the  case  of  the  spline  ADI  method,  the  boundary  conditions 
have  to  be  imposed  in  such  a  way  that,  at  convergence,  they  are  to  be  fourth- 
order-accurate,  consistently  with  the  numerical  scheme.  Dirichlet  boundary 
conditions  are  obvious,  insofar  as,  if  to  or  \J;  are  prescribed  at  the  boundary, 
it  is  required  that 

Di)/  =  Dw»D4>  =  Du)  =  0  (17) 

at  the  appropriate  boundary  gridpoints.  A  Neumann  boundary  condition  for  the 
stream  function  is  slightly  more  complicated  to  deal  with.  Here  only  one 
example  will  be  given,  namely  how  to  impose  the  boundary  condition  ^  =  0 


in  the  first  sweep  of  the  spline  ADI  method:  The  boundary  condition 
itself  is  written  in  incremental  form  and  replaced  by  its  fourth-order- 
accurate  spline  approximation  (Eqn.  4);  the  stream  function  equation  at 
the  boundary  and  the  second  boundary  condition  on  the  stream  funct i on  are 
also  considered.  In  this  way,  three  linear  equations  are  available  for 
evaluating  the  vorticity  and  the  stream  function  at  the  boundary  point  and 
the  stream  function  at  a  mirror-image  gridpoint (external  to  the  flow  field). 
The  stream  function  equation  basically  allows  to  eliminate  this  extra  un¬ 
known  and  the  two  boundary  conditions  on  the  stream  function  allow  a  direct 
evaluation  of  the  wall  stream  function  and  vorticity.  Obviously,  all  the 
K  terms  are  known  from  the  previous  time  level  and  the  value  corresponding 
to  the  aforementioned  mirror-image  gridpoint  is  obtained  by  linear  or 
quadratic  extrapolation.  Several  other  boundary  conditions  are  possible; 
for  example,  using  alternate  expressions  of  m  and  M,  which  do  not  introduce 
a  mirror- image  gridpoint,  or  two  different  expressions  for  m,  which  amounts 
to  enforce  that  the  spline  interpolating  polynomial  has  a  continuous  first 
derivative  thru  the  boundary  gridpoint.  In  the  present  study  all  the 
approaches  described  above  have  been  used  successfully. 

In  Ref.  1,  an  alternate,  Crank  Nicolson-type,  linearization  and  dis¬ 
cretization  in  time  of  the  governing  equations  was  also  used.  The  present 
techniques  also  have  this  option:  In  the  computer  programs  all  coefficients 
of  the  block-tridiagonal  systems,  to  be  solved  at  every  sweep  of  the  ADI 
methods,  contain  a  CR  coefficient  which  can  be  either  1  or  0.5,  to  provide 
the  implicit  backward  (see,  e.g.,  Eqns.  9  and  10)  or  the  Crank  Nicolson 
time  discretization,  directly.  The  Crank  Nicolson  approach  is  obviously 
second-order-accurate  in  time,  but,  since  the  present  techniques  only  provide 
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the  steady  state  solution,  this  is  not  necessarily  an  advantage.  Further¬ 
more,  when  using  a  Crank  Nicolson  averaging,  in  order  to  obtain  a  correct 
value  of  the  vorticlty  at  the  wall,  the  boundary  conditions  have  to  be 
inconsistent  with  the  difference  equations;  that  is,  it  is  necessary  to  use 
an  implicit  backward  time  discretization  of  the  stream  function  equation  at 
all  boundary  points  in  the  ^  sweep  of  the  ADI  methods,  in  order  to  obtain  the 

correct  vorticity  at  the  wall  at  convergence.  This  result  is  consistent  with 

9  10 

that  already  observed  by  Davis  et  al  and  Briley  and  McDonald  .  For  this 
reason,  most  of  the  results  later  presented  in  this  report  have  been  obtained 
with  the  fully  implicit  time  discretization  of  Eqns.  (9,10).  However,  in  some 
calculations,  the  Crank  Nicolson  time  linearization  (which  obviously  produces 
identical  results  at  convergence)  has  been  found  to  provide  faster  convergence 
rates  for  the  same  values  of  At. 
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SECTION  IV 


RESULTS 


A.  Flow  Between  Two  Rotating  Circles 

The  present  methods  have  been  developed  in  order  to  compute  viscous 
steady  flows  past  two-dimensional  airfoils,  in  connection  with  a  method 
for  generating  a  system  or  orthogonal  curvilinear  coordinates.  The  viscous 
flow  between  two  concentric  rotating  circles  has  been  considered  as  a  model 
problem  (somewhat  simulating  such  a  flow  configuration),  for  which  an  exact 
solution  is  available  for  comparisons.  The  inner  circle  of  radius  equal  to 
one  is  chosen  to  be  stationary  in  order  to  test  the  no-slip,  zero  injection 
boundary  conditions,  usually  given  at  the  surface  of  a  stationary  airfoil, 
whereas  the  outer  circle  rotates  at  such  a  speed  that  the  vorticity  on  its 
boundary  is  also  equal  to  one.  A  given  vorticity  has  been  imposed  at  the 
outer  circle  (of  radius  equal  two),  because,  for  external  flow  configurations, 
the  outer  boundary  is  usually  chosen  at  a  sufficient  distance  from  the 
surface  of  the  airfoil,  that  a  zero  vorticity  boundary  condition  is  imposed. 
The  physical  flow  field,  divided  into  a  system  of  equally  spaced  polar 
coordinates,  has  been  transformed  into  a  rectangle  in  the  £,r|  plane.  All 
metric  coefficients  and  the  Jacobian  have  then  been  evaluated  numerically 
with  fourth-order-accurate  spline  interpolating  polynomials,  except  for  the 
mixed  derivatives  which  are  identically  zero.  It  is  noteworthy  that,  due 
to  the  coordinate  transformation,  the  flow  field  is  "opened  up,"  so  that 
periodic  boundary  conditions  are  needed  in  the  £  direction.  For  more  details 
and  for  an  algorithm  capable  of  solving  periodic  2x2  block-tridiagonal 
systems,  the  reader  is  referred  to  Ref.  1. 
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The  results  obtained  with  both  the  ADI  and  the  spline  ADI  methods  for 
such  a  test  problem  are  given  in  Fig.  1,  where  the  stream  function  at  the 
center  of  the  anulus  and  the  vorticity  at  the  wall  of  the  inner  circle  are 
plotted  versus  the  inverse  of  the  number  of  gridpoints  (A),  square  (for  the 
ADI  method),  and  to  the  fourth  power  (for  the  spline  ADI  method),  for  several 
values  of  A.  Two  sets  of  spline  ADI  results  are  given,  corresponding  to  the 
use  of  a  linear  and  a  quadratic  extrapolations  for  the  K  spline  functions. 

All  results  are  seen  to  tend  to  the  exact  solution  as  A  tends  to  zero  and 
with  the  correct  second-order-accuracy  and  fourth-order-accuracy,  respectively. 
Also,  the  higher  order  extrapolation  is  seen  to  produce  more  accurate  results 
as  expected.  However,  the  most  interesting  point  to  make  is  that,  whereas  the 
computation  cost  is  almost  equivalent  (the  convergence  rate  is  almost  the 
same  for  both  approaches,  probably  due  to  the  one-dimensional  nature  of  the 
problem),  the  more  accurate  spline  ADI  results,  using  a  10  x  10  mesh,  have  a 
discretization  error  up  to  five  times  smaller  than  that  of  the  ADI  results, 
using  a  24  x  24  mesh. 

B.  The  Driven  Cavity  Flow 

The  classical  driven  cavity  problem  (see,  e.g.,  Rubin”'"'*’)  was  also  used  as 
a  test  problem  to  verify  the  present  numerical  techniques.  The  Re  =  100  case 
has  been  considered,  using  a  uniform  rather  coarse  14  x  14  mesh.  The  values 
of  the  vorticity  at  the  center  of  the  moving  wall  of  the  cavity  and  the  maximum 
value  of  the  stream  function  are  given  in  Table  1  for  the  ADI  method  as  well  as 
for  the  spline  ADI  method  using  either  linear  or  quadratic  extrapolation  for 
the  spline  boundary  conditions.  The  second-order-accurate  results  are  seen  to 
perfectly  coincide  with  the  results  given  by  Rubin*"*".  It  is  worth  mentioning 
that  the  present  approach  does  not  use  the  divergence  form  of  the  equations. 


The  results  in  Table  1  have  been  obtained  with  both  the  backward  and  Crank 


Nicolson  time  discretization  of  the  governing  equations  and  with  values  of 
A  and  B  equal  to  or  greater  than  one.  It  is  noteworthy  that,  if  one  uses 
CR  =  0.5  and  the  consistent  boundary  conditions,  the  vorticity  at  the  wall 
is  completely  wrong  although  a  correct  solution  is  obtained  at  all  internal 
gridpoints.  The  spline  ADI  results  given  in  Table  1  clearly  indicate  the 
superior  accuracy  obtained  by  means  of  the  high  order  spline  correction 
procedure  (see  Ref.  11  for  very  accurate  results).  It  is  noteworthy  that 
in  the  present  calculations  it  was  unnecessary  to  use  values  of  A  and  B 
greater  than  1.  The  results  in  Table  1  have  been  obtained  within  2  CPU 
minutes  on  an  HP  1000/F  minicomputer,  thus  verifying  the  efficiency  of  the 
proposed  methods.  The  Re  =  1000  rather  difficult  case  has  also  been  con¬ 
sidered,  using  a  uniform  20  x  20  mesh,  and  the  results  obtained  with  the  ADI 
method  are  given  in  Table  1.  Although  insufficiently  accurate,  due  to  the 
use  of  a  uniform  mesh,  which  is  completely  incapable  of  capturing  the  thin 
boundary  layers  at  the  walls  of  the  cavity,  the  results  in  Table  1  are  the 
correct  ones  corresponding  to  a  second-order-accurate  central  space  discreti¬ 
zation  of  the  steady  state  equations.  Furthermore,  the  most  important  point 
is  that,  whereas  this  Re  =  1000  case  was  found  to  be  an  impossible  task  for 
the  method  of  Ref.  1,  a  fully  converged  solution  was  obtained  in  less  than 
1000  ADI  cycles  (30  CPU  minutes  on  the  HP/1000F  minicomputer)  by  using  the 
present  improved  method. 

C.  Flow  in  a  Channel  of  Complex  Geometry 

The  present  methods  were  finally  used  to  compute  viscous  laminar  flow 

inside  a  channel  of  complex  geometry.  The  problem  was  first  proposed  by 
12 

Roache  who  numerically  verified  that,  if  the  length  of  the  channel  is 
scaled  proportionally  to  the  Reynolds  number  of  the  flow,  self  similar  flow 
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conditions  are  obtained  for  very  high  Re  values.  Recently  the  same  problem 
has  been  used  as  a  numerical  test-case  for  comparing  the  accuracy  and  effi¬ 
ciency  of  several  numerical  Navier-Stokes  solvers  by  the  IAHR  working  group 
on  refined  modelling  of  flows  in  its  VI  meeting  held  in  Rome  (June  24-25, 
1982).  The  geometry  of  the  channel  is  given  in  Figs.  2a  for  the  Re  =  10  case 
and  2b  for  the  Re  =  100  case.  The  lower  wall  of  the  channel  is  given  ana¬ 


lytically  as 


Y^  =  -|  [tanh  (2  -  30  x  /Re)  -  tanh  2] 


and  its  centerline  as 


Y  =  1 
u 


The  inlet  and  outlet  sections  of  the  channel  are  finally  given  as 


x  =  0  and  x  =  Re/ 3 


(20a, b) 


respectively. 

In  the  present  study,  a  system  of  orthogonal,  curvilinear  coordinates 

has  been  used  to  map  the  physical  (x,y)  flow  domain  into  a  rectangle  in  the 

(C,n)  computational  domain.  A  simple  algebraic  transformation,  as  given  by 
13  14 

Blottner  and  Ellis  and  described  by  Davis  ,  has  been  used  for  simplicity 
as  well  as  for  taking  full  advantage  of  the  shape  of  the  channel,  being  pre¬ 
scribed  analytically.  The  system  of  (r|  =  constant)  coordinate  lines  in  the 
physical  plane  has  been  prescribed  as  follows:  The  line  n  =  0  coincides  with 
the  lower  boundary  (the  wall)  of  the  channel  and  the  line  n  =  1  with  its 
upper  boundary  (the  symmetry  line).  All  the  other  r)  =  (j  -  l)Ar|  (j  =  2,3,... 
N+l,  An  *  ^)  are  given  by  the  following  expressions 


N  j-1  j-1 

— 1_3 -  y  (x)  +  ^ - J— i  Y  (x) 

N  .  Vx>  N  ,  uw 
q  -  1  q  -  1 


I )  i  ■* 


A  family  of  (5  m  constant)  lines,  orthogonal  to  the  r)  =  constant  lines  have 
then  been  obtained  by  integrating  numerically  with  a  fourth-order-accurate 
Runge  Kutta  procedure,  the  following  equation 

rh.\  (iZ) 

i*  _  x  .... 

3n  "  -  $>; 

starting  from  prescribed  points  at  the  lower  boundary  of  the  channel.  A  few 
points  are  of  interest:  The  distance  between  two  successive  0  =  constant 


lines  in  the  physical  plane  has  been  chosen  to  increase  at  a  constant  rate 
(q  =  1.1)  starting  from  the  lower  boundary.  In  this  way  a  finer  resolution 
is  obtained  near  the  wall  of  the  channel  where  viscous  effects  are  more 
important.  The  distance  between  two  successive  5  =  constant  lines  along 
the  x  coordinate  in  the  physical  plane  has  also  been  chosen  to  increase  at  a 
constant  rate  (r  =  1.054),  starting  from  the  point  x^-  Re/ 15  (where  the  lower 
boundary  of  the  channel  has  an  inflection  point)  in  both  x  >  x^  and  x  <  x^ 
directions.  In  this  way  a  finer  resolution  is  obtained  in  the  region  where 
a  separation  bubble  is  likely  to  develop.  Finally,  the  5  =  0  and  5=1  lines 
in  the  physical  plane  have  been  chosen  to  coincide  with  the  entrance  and  the 
exit  of  the  channel,  i.e.,  with  the  x  =  0  and  x  *  Re/3  lines,  for  convenience, 
and  are  not  perfectly  orthogonal  to  the  p  =  constant  lines.  Therefore,  the 
metric  coefficients  multiplying  the  mixed  derivatives  in  the  governing  Navier- 
Stokes  equations  are  not  identically  zero  at  all  gridpoints  and  the  accuracy 
of  the  spline  ADI  method  slightly  deteriorates  locally,  insofar  as  the  mixed 
derivatives  are  evaluated  with  standard  second-order-accurate  central  diff¬ 
erences.  The  curvilinear  coordinates,  generated  as  described  above,  are 


shown  in  Fig.  3  for  the  Re  *  10  case,  together  with  the  boundary  conditions 
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used  in  the  present  calculations.  After  evaluating  all  the  grldpoint  locations 
in  the  physical  x,y  plane  (corresponding  to  a  uniform  cartesian  grid  in  the 
computational  £,n  plane),  the  metric  coefficients  a,  B,  Y,  a,  T,  and  the  jacobian 
of  the  transformation  are  evaluated  at  all  internal  gridpoints  by  means  of 
central  differences  or  fourth-order-accurate  spline  approximations  in  the  ADI 
and  spline  ADI  methods,  respectively.  At  the  boundary  points  the  metric  co¬ 
efficients  necessary  to  evaluate  oo  from  the  stream  function  equation  are 
obtained  by  linear  or  cubic  extrapolation  from  the  neighboring  gridpoints. 

The  numerical  solutions  obtained  for  the  Re  =  10  and  the  Re  =  100  cases  by  means 
of  both  present  ADI  methods  are  given  in  Table  2  as  the  values  of  the  vorticity 

at  the  wall  versus  the  x.  locations  of  the  L  coordinate  lines  (the  x.  values 

i  l 

corresponding  to  Re  =  100  have  to  be  multiplied  by  10).  For  convenience,  these 
results  are  also  plotted  in  Figs.  4  and  5  for  the  Re  =  10  and  Re  =  100  cases, 
respectively;  the  results  obtained  by  means  of  the  spline  ADI  method,  using  a 
coarser  10  x  10  mesh,  are  also  given.  All  solutions  are  seen  to  coincide,  for 
all  practical  purposes,  and  favorably  compare  with  those  obtained  by  means  of 
several  other  methods  (see  e.g..  Ref.  15).  The  20  x  20  mesh  spline  ADI  solution 
is  the  most  accurate  and  the  other  two  solutions  have  comparable  accuracy.  For 
the  two  ADI  calculations,  the  minimum  value  of  the  vorticity  at  the  wall  is 
plotted  versus  the  normalized  number  of  iterations,  to  provide  an  idea  of  the 
convergence  properties  of  the  technique.  A  full  convergence  (to  machine  accu¬ 
racy)  has  been  obtained  within  about  130  and  150  iterations  (At  =  0.075  and 
At  =  0.08)  for  the  Re  =  10  and  Re  *  100  cases,  corresponding  to  less  than  5  CPU 
minutes  on  the  HP  1000/F  minicomputer.  If  one  considers  that  a  reasonable  con¬ 
vergence  is  obtained  in  about  40%  of  the  total  number  of  iterations  (see 
Figs.  6  and  7),  the  efficiency  of  the  present  ADI  approach  is  self  evident. 

The  spline  ADI  procedure,  using  the  same  mesh,  required  from  five  to  ten  tiroes 
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more  iterations  to  fully  converge,  because  a  stability  limitation  on  At  had 


to  be  satisfied  for  A  =  B  =  1.  An  optimization  of  the  convergence,  obtained 
by  using  different  values  of  A  and  B,  or  by  correcting  the  ADI  solution  only 
after  a  given  number  of  ADI  cycles,  has  not  been  pursued.  However,  by  com¬ 
paring  the  ADI  solutions  with  the  10  x  10  mesh  spline  ADI  solutions  (which 
have  comparable  accuracy) ,  it  turns  out  that  these  required  about  20%  less 
CPU  time. 
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SECTION  V 


CONCLUSIONS 

An  improved  ADI  numerical  technique  for  the  solution  of  incompressible 
viscous  steady  flows  has  been  developed,  together  with  a  fourth-order-accurate 
spline  ADI  method  obtained  by  applying  a  spline  deferred  corrector  approach 
to  the  present  ADI  technique.  The  validity  and  efficiency  of  the  present 
approaches  have  been  demonstrated  by  their  application  to  the  numerical 
solution  of  three  viscous  flow  problems. 
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SADI  din. 


TABLE  1 


DRIVEN  CAVITY  RESULTS 


Re  «  100  14  x  14  Mesh 


‘Vtp 

*MAX 

ADI 

-8.196 

-0.0874 

Spline  ADI 

-8.557 

-0.0941 

(linear  extrap.) 

Spline  ADI 

/ _ l  — 

-9.380 

-0.0992 

(quad,  extrap.) 


Re  =  1000 

20  x  20  Mesh 

wMrP 

^MAX 

ADI 

-28.471 

-0.0344 
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y  =  O,  dty/Zn  =  o 

Figure  3.  Re  =  lO  Computational  Grid  and  Boundary  Conditions 


TABLE  2 


! 


CHANNEL  FLOW  ADI  AND  SPLINE  ADI  WALL  VORTICITY  RESULTS 


i 

X 

1 

0.0 

2 

0.1477 

3 

0.2879 

4 

0.4209 

5 

0.5470 

6 

0.6666 

7 

0.7862 

8 

0.9123 

9 

1.0453 

10 

1.1855 

11 

1.3333 

12 

1.4891 

13 

1.6534 

14 

1.8266 

15 

2.0092 

16 

2.2018 

17 

2.4048 

18 

2.6188 

19 

2.8455 

20 

3.0824 

21 

-  3.3333 

SADI 


Re*5 10 

Re=100 

3.0756 

3.0769 

2.5877 

2.5393 

2.1025 

1.9546 

0.9908 

1.1806 

0.1826 

0.4915 

-0.0864 

0.0992 

-0.1314 

-0.0679 

-0.1145 

-0.1227 

-0.1024 

-0.1261 

-0.1025 

-0.0930 

-0.0986 

-0.0349 

-0.0752 

0.0431 

-0.0292 

0.1281 

0.0337 

0.2158 

0.1056 

0.2966 

0.1798 

0.3747 

0.2511 

0.4399 

0.3147 

0.5020 

0.3661 

0.5507 

0.4001 

0.5945 

0.4115 

0.6415 

ADI 


Re=10 

Re=100 

3.0113 

2.9970 

2.6612 

2.5620 

2.0040 

1.9564 

0.9675 

1.1820 

0.2148 

0.4999 

-0.0837 

0.1075 

-0.1365 

-0.0621 

-0.1184 

-0.1189 

-0.1037 

-0.1245 

-0.1032 

-0.0901 

-0.0993 

-0.0327 

-0.0761 

0.0482 

-0.0307 

0.1307 

0.0313 

0.2212 

0.1024 

0.2985 

0.1762 

0.3790 

0.2473 

0.4396 

0.3114 

0.5058 

0.3636 

0.5470 

0.3986 

0.5997 

0.4336 

0.6525 

Figure  4.  Re  =  10  Wall  Vortlcity 
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Figure  5.  Re  =  100  Wall  Vorticity 


Figure  6.  Minimum  Wall  Vorticity  Versus  Normalized  Iteration  Number  (ADI,  Re  -  lO) 
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